library(Seurat)
library(tidyseurat)
library(harmony)
library(SingleR)
library(ggpubr)
library(tidyverse)
source('00_util_scripts/mod_seurat.R')

fq.come <-
  list.files('/home/supervisor/mist2/gjsx/ye2024_AT_cOME/',
             full.names = T, recursive = T, pattern = '\\d.fq.gz')

tibble(fq = fq.come) |> 
  mutate(sample = str_extract(fq, '.{5}(?=_r)'),
         pair = ifelse(str_detect(fq, 'r1'), 'fastq_1', 'fastq_2')) |>
  pivot_wider(names_from = 'pair', values_from = 'fq') |> 
  unnest(-1) -> nf.come

nf.come |>
  mutate(readgroup = sample) |>
  relocate(sample, readgroup) |>
  write_csv('~/append-ssd/nextflowing/ye.ah.come1/binext.csv')

# fastp vs raw
fastped <- Read10X('~/append-ssd/nextflowing/ye.ah.come1/M1025.fastpSolo.out/Gene/filtered/')

sobj_fp <- fastped |>
  CreateSeuratObject(min.cells = 3, min.features = 200)

sobj_fp <- sobj_fp |>
  PercentageFeatureSet('^MT-', col.name = 'mito.ratio')

g1 <- sobj_fp |> ggplot(aes(mito.ratio)) +
  geom_boxplot(color = 'blue') +
  expand_limits(x = 31)

raw <- Read10X('~/append-ssd/nextflowing/ye.ah.come1/M1025Solo.out/Gene/filtered/')

sobj_rw <- raw |>
  CreateSeuratObject(min.cells = 3, min.features = 200) |>
  PercentageFeatureSet('^MT-', col.name = 'mito.ratio')

g2 <- sobj_rw |> ggplot(aes(mito.ratio)) +
  geom_boxplot(color = 'red') +
  expand_limits(x = 31)

g1 / g2

sobj_rw <- sobj_rw |> quick_process_seurat()

sobj_fp <- quick_process_seurat(sobj_fp)

# all star -------
all.star <- list.files('~/append-ssd/nextflowing/ye.ah.come1/empty_intron/',
                       full.names = T,
           include.dirs = T, recursive = T, pattern = 'filtered') |>
  str_subset('filtered_full')

ah9 <- all.star |>
  str_extract('[A-Z].{4}(?=/Solo)') |>
  set_names(all.star, nm = _) |>
  Read10X()

sobj <- ah9 |>
  CreateSeuratObject(min.cells = 3, min.features = 200)

sobj <- sobj |>
  PercentageFeatureSet('^MT-', col.name = 'mito.ratio')

sobj |> VlnPlot('mito.ratio',pt.size = 0)

sobj <- sobj |>
  filter(mito.ratio < 10)

sobj <- sobj |>
  quick_process_seurat(batch = c('orig.ident', 'sex'), leiden = F)

# infer sample sex from XIST expr
sobj |>
  DotPlot(c('XIST'), group.by = 'orig.ident')

ah.meta <- read_delim('mission/ye_AH-cOME/ah-come.meta.txt')

# compute odds ratio --------
ah.meta |>
  mutate(i2t = ifelse(`FCGR2B-I232T` == 'II', 0, 1),
         g2r = ifelse(`IGHG1-G396R` == 'GG', 0, 1),
         sex = ifelse(sex == 'male', 0, 1)) |>
  glm(cOME ~ i2t + sex, family = 'binomial', data = _) |>
  broom::tidy()

ah.meta |>
  ggplot(aes(`FCGR2B-I232T`, fill = as_factor(cOME))) +
  geom_bar() +
  theme_pubr(legend = 'none') +
  scale_fill_manual(values = c('blue','red'),
                    labels = c('AH','AH + chronic OME')) +
  labs(fill = 'Group')

g1 <- last_plot()

ah.meta |>
  ggplot(aes(`IGHG1-G396R`, fill = as_factor(cOME))) +
  geom_bar() +
  theme_pubr(legend = 'right') +
  scale_fill_manual(values = c('blue','red'),
                    labels = c('AH (n=5)','AH + chronic OME (n=4)')) +
  labs(fill = 'Group')

g2 <- last_plot()

g1 + g2 + patchwork::plot_layout(widths = c(2:3))

sobj <- ah.meta |>
  mutate(orig.ident = id) |>
  left_join(sobj, y = _)

sobj |> write_rds('mission/ye_AH-cOME/ah9.rds')

hpca.imm <- read_rds('00_util_scripts/ref/hpca_immune.rds')

sobj <- sobj |>
  mark_cell_type_singler(ref = hpca.imm,
                         fine_label = T,
                         new_label = 'hpca.fine')

sobj |>
  DimPlot(group.by = 'hpca.fine',
          cols = 'Paired')

sobj |>
  DotPlot(seurat_markers) +
  RotatedAxis()

